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: Abstract 

We give a rigorous derivation of Fourier's law from a system of closure equations for a 
nonequilibrium stationary state of a Hamiltonian system of coupled oscillators subjected 
to heat baths on the boundary. The local heat flux is proportional to the temperature 
^ I gradient with a temperature dependent heat conductivity and the stationary temperature 

■ exhibits a nonlinear profile. 

o 

a^ 
o 
\o 

^ i One of the simplest and most fundamental nonequilibrium phenomena is heat conduction 

■ in solids. It is described by a macroscopic equation, Fourier's law, which states that a local 
I . temperature gradient is associated with a flux of heat which is proportional to the gradient: 

: J{x) = -k{T{x))VT{x 



where the heat conductivity k{T{x)) is a function only of the temperature at x. 

Despite its fundamental nature, a derivation of Fourier's law from first principles, or even 
^ ■ within a suitable approximation, such as a Boltzmann type equation, lies well beyond what 

can be mathematically proven (for reviews on the status of this problem, see jS], [HI and [13j). 
The quantities T and JT" in are macroscopic variables, statistical averages of the variables 
describing the microscopic dynamics of matter. A first principle derivation of entails a 
definition of T and in terms of the microscopic variables and a proof of the law in some 
appropriate limit. 

In this letter we outline a rigorous proof ,4J of Fourier's law starting from a closure ap- 
proximation of the equations for the nonequilibrium stationary state of a Hamiltonian system 
subjected to boundary heat baths. The physical situation we have in mind is a slab of crystal 
of linear extension heated at the two ends by temperatures Ti and T2. In this case one would 
expect VT and J' to be 0(1/A^) and ((T)) to hold up to corrections of order o{l/N). This is 
indeed what we prove in our model together with a detailed description of the temperature 
distribution in the bulk. 

Since dH) is a macroscopic law it is expected to hold for a classical system as well as for a 
quantum one and the quantum corrections are expected to be small except at low temperatures. 
A classical toy model describing the above situation which has been intensively discussed in 
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recent years is given by coupled oscillators organized on a strip of width N in d-dimensional 
cubic lattice Z'^. The oscillators are indexed by lattice points x = (xi, . . . , Xd) with < xi < 
and carry momenta and coordinates {pxiQx)- The dynamics of the oscillators consists of two 
parts: Hamiltonian dynamics in the bulk and noise on the boundary modelling heat baths at 
temperatures Ti and T2. 

We consider a Hamiltonian of the form 

H{q,p) = \Y.Pl + \i<l^^"<l) + (2) 

which describes a system of coupled anharmonic oscillators with coupling matrix a;^, i.e. 
(g, cu^g) = J^QxQy^'^ix — y). The noise is specified by the Gaussian random variables ^xit) 
at sites x on the boundary with covariance 

< Ut)Ut') >= '^iS.yiTidx.o + T26x,N)6{t - t') = 2Cxy6{t - t'). (3) 
The (stochastic) dynamics is given by (jx = Px and 

Px = i-Q^ - IxPx) + ix (4) 

where the friction is 7^: = 7(^x10 + ^^lAf) (more precisely, (jH) is a Ito stochastic differential 
equation). These equations define a Markov process {q{t),p{t)) and we are interested in the 
stationary states for this process. 

In the equilibrium case of equal temperatures, Ti = T2 = T, an explicit stationary state 
is given by the Gibbs state Z~^e~^^^'''^^dqdp of the Hamiltonian H with inverse temperature 
/? = 1/T. When Ti 7^ T2 there is no such simple formula and, indeed, in our setup, even the 
existence of a stationary state is an open problem. In the d = 1 case of a finite chain of 
oscillators the existence is proved under conditions (see for a review of such results, first 
obtained in jH]) as well as the convergence of the Markov process to this state as t — > 00. 

Supposing that we have a stationary state, let us formulate the statement (^. Writing H 
as a sum of local terms, each one pertaining to a single oscillator: H = J^xeA^x, one has up 
to noise terms = V ■ j(x), where the microscopic heat current j{x) will depend on py and 
qy for y near x provided the coupling matrix cu^ is short ranged. Let also t{x) = ^pl be the 
kinetic energy of the oscillator indexed by x. Then, the macroscopic temperature and heat 
current in eq. (P) are defined by T{x) =< t{x) >^ and J{x) =< j{x) >^ where < ■ >^ denotes 
expectation in the stationary state. 

The only rigorous results in our model as far as deriving are for the harmonic case 
of quadratic H ^TJ ^]. In that case, Fourier's law does not hold: the current j{x) is 0{1) 
as N ^ 00 whereas VT = except near the boundary. If A 7^ the law seems to hold in 
simulations in all dimensions pp. In an analogous momentum conserving model conductivity 
seems anomalous in low dimensions: in depends on as A^" in d = 1 and logarithmically 
in d = 2. It is a major challenge to explain the a which, numerically, seems to be in the interval 
[1/3,2/5] (see [Zj, 0, ^ for theories on a). 

In this paper we suppose the stationary state exists and we study its properties via its 
correlation functions. Let us denote {qx,Px) = {uix,U2x), A(u®^)q,^. = —X6a,2(ll, (Xu)x = 
{0,'jxPx)^ and 7] = (0,^)^. Then (jH) becomes 

u{t) = [{A-r)u + A{u'^^))+r]{t) (5) 
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where A = ^2 g j • Applying (0) to the stationary state correlation functions 

Gn{Xi, ...,Xn) =< Ux'i ® ... ® > 

these are seen to be solutions of the Hopf equations 

(A„ - r„)G„ + A„G„+2 + C„G„_2 = 0. (6) 

where An = J2 Axi and r„ is defined similarly. A„ and C„ are linear operators involving the 
quartic vertex and noise covariance Cxy of respectively. The equations (jHl) have the drawback 
that they do not "close": to solve for G„, we need to know Gn+2- 

We will now introduce an approximation that will lead to a closed set of nonlinear equations 
for G2- For this we note that for A small, the equilibrium Ti = T2 measure is close to Gaussian. 
When Ti 7^ T2 we expect this to remain true and we look for a Gaussian approximation to 
equation ® for small A by means of a closure, i.e. expressing the G„ in terms of G2. Let 
G4 be the connected correlation function describing deviation from Gaussianity. Then the first 
equation in the hierarchy ^ reads: 

{A2 -T2 + S2)G2 + A2GI + C = (7) 

where 112(^2)^2 = A2J2G2 ® G2. The simplest closure would be to drop the A2GI from 0. 
This leads to a nonlinear equation for G2. It turns out that the solution to this equation is 
qualitatively similar to the A = case, i.e. G2 does not exhibit a temperature profile nor a 
finite conductivity. The only effect of the nonlinearity is a renormalization of uj. 
The next equation in the hierarchy becomes after some some algebra 

(A4 - r4 + S4)G^ + b{G2) + AiGl = 0, (8) 

where Gg is the connected six point function, 

^4{G2)Gl = Y.iMG2 ® Gl) -G2® A2G'2). 

and 

6(G2) =EA1(G'2®G2®G2), 

p 

where A4 has A acting on all the three factors 6*2. 

((Zj) and (jHI) yield an exact equation for the two point function with the connected six point 
function as an input. Our closure approximation consists of dropping the Gg term in (jH)) thereby 
yielding a closed set of equations for G2. For simplicity we also drop the operators S2, S4 and 
r4: these could be included in our analysis, but do not change the main structure that is due 
to the term h{G2)- Hence the closure equation we study is for G = G2: 

{A2~T2)G + Af{G)+C = Q (9) 

with 

AT = -A2 A^'b{G). (10) 
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we can write G 



To write this more concretely, let us introduce the matrices Qxy =< QxQy >, Pxy =< PxPy > 
and Jxy =< QxPy >■ Clearly Q = J + so the (l,l)-component of says Jxy = —Jyx and 

Q J 
-J P 

Next, anticipating translation invariance in the directions orthogonal to the 1-direction we 
write 

G{x,y) = f e'P'^'''+y'^+'''^''-y^G{p,k)dpdk (11) 



where the integrals over p and ki are Riemann sums on a 2^ lattice. The inverse of ^4 is 
written as 

-A-'= / e*^Mt = / mf^dt 
Jo Jo 

where R{t) = e*^. In Fourier space the latter is 

), (12) 



2 ^f^^ V isuj{(l) 

Then some algebra yields the following expression for the nonlinear term 
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Nip, k) = J2j MT. ^i^iPi + ^i) + ' n WsAPi, ki)-(^l is4uj{p4 + ki) 

ssooips + ks) [uj{p3 + k3y^6{2p3)Ws,{p4, k^) - uj{pi + ki)-H{2pi)Ws,{pz, A;3)] (13) 

where 

du = 6{2p - J2{p, + h))6(£{p, - h))6{p -k~p^- k^)dpdk, (14) 
and k = (A;i)f=i and similarily for p and s. W is the following combination 

Ws{p,q) = Qip,q) + isuj{p + qy^J{p,q). (15) 

Eq. Q is a nonlinear set of equations for the pair correlation functions of our model. In 
|3] we have proven that they have a unique solution which we now proceed to explain. The 
1,2-component of Q (coming from ^ < qp >) gives 

P = uj{p, kfQ + i(( jr - r J) - Ni2{p, k) - Ni2{p, -k)) (16) 

where Lj{p,k)^ = ^{uj{p + kY + uj{p — kY). Since the nonlinear term M depends only on Q 
and J this expresses P in terms of them. The rest of then yields two equations for the two 
unknown functions Q and J which become 

Su^iQ, jf + ^^iQ, J) + (r J + jr, tp + pvf = (o, c)^ (i7) 

where N{Q, J) = {Af^ip, -k) - Mi2{p, k), -Af22ip, k)Y and 5uo'^{p, k) = uj{p + kf - uo{p - kf . 
An important property of M is that, for all T, M{TQq, 0) = where Qo = For 7 = 
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these form a 1-parameter family of solutions of (fTTj) . In the equilibrium case Ti = T2 and 7 7^ 
only one of these persists, namely the one with T = Ti. This is the analogue in the closure 
equation of the true equilibrium Gibbs state which has Q = Qo + 0{X). 

Since we are looking for a solution that is locally in xi close to this equilibrium we need to 
understand the linearization of the nonlinear term A/" at Go- It turns out that it is given by an 
operator which is a multiplication operator in the slow variable p: 

^{TQo + 6Q, 6J){p, k) = C,i6Jip, ■), 6Q{p, ■)f{k). (18) 

£p is a matrix of operators Cij{p). Each of these acts on functions of as a sum of a multipli- 
cation and an integral operator 

{C,,{p)f){k) = A{p, k)f{k) + J B{p, k, k')f{k')dk'. (19) 

The integral kernel B{p, k, k') is of the form 

B{p, k,k') = = ^M^i) + ^Mk' + P) + SiUj{p - k)) 

s i=l 

■5{k — ki — k2 — k')ps{ki, k2, k', k,p)dkidk2 (20) 

where A(x) = 6{x) or V (^^^ and ps is a smooth function. A{p, k) is given by a similar expression 
integrated over k'. 

The integrand in (f^UI) represents phonon scattering. The delta functions impose momentum 
and energy conservation when p = (it turns out that only terms with Y^Si = contribute) 
which is the point where our functions are peaked: the translation invariant equilibrium has 
support at p = and the nonequilibrium solution will also have most of its mass in the 
neighborhood of this point. Thus it is important to understand Co- For parity reasons Cij{0) = 
for i 7^ j whereas i2ii(0) is invertible. Invertibility of Cq would then follow from invertibility 
of £22(0). This, however, is not the case: £22(0) has two zero modes. 

One of them is easy to understand. Since M{TQq, 0) = for all T, taking derivative 
with respect to T, one finds £22(0)0;^^ = 0. There is, however, also a second zero mode: 
£22(0)cj-3 _ 0. 

While the first zero mode has to persist for the full Hopf equations due to the one parameter 
family of Gibbs states that solve them for 7 = 0, the second one is an artifact of the closure 
approximation. The phonon scattering described by the nonlinear term conserves phonon 
energy, leading to the first zero mode, and also phonon number, leading to the second one. 
The connected six-point correlation function which was neglected in the closure approximation 
would produce terms that violate phonon number conservation and remove the second zero 
mode. However, for weak anharmonicity its eigenvalue would be close to zero and should be 
treated as some perturbation of the present analysis. 

The second zero mode leads one to expect that our equations have in the 7 = limit a two 
parameter family of stationary solutions which indeed is the case. These are given by 

QtA^.v) =T J e'^^''-y\uj{kf - Auj{k))-^dk. (21) 

The second zero-mode is proportional to the derivative of Qt,a with respect to A, at A = 0. 
We are then led to look for solutions in the form 

Q(x, y) = Qt{^),a{^){x -y) + q{x, y), (22) 
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where the first term is of local equilibrium form with slowly varying temperature and "chemical 
potential" profiles T(x) and ^(x) and where g is a perturbation orthogonal to the zero modes 
in a suitable inner product. 

Projecting equation (|T7jl on the complementary subspace of the zero modes yields a nonlin- 
ear equation for J and q with an invertible linear part. It can be solved by fixed point methods 
in a suitable Banach space and yields J and q as functionals of T and A. The heat and phonon 
number currents are given in momentum space by 



for a = 1, respectively and are thus nonlinear functionals of VT and VA. This relation is 
the precise form of the Fourier law. In particular we get relation (0) for the heat current, up 

to corrections 0{^^^^) (here AT = T2 — Ti), with the thermal conductivity given by: 



Projecting then equation p7|) to the two left zero modes of £22(0) yields two conservation 
laws for the currents which become nonlinear (and nonlocal) elliptic equations for the functions 
T(x) and A{x) whose solutions give for the inverse temperature /3(x) = T{x)~^ a linear profile: 

P{x) = Pi + jj-{P2 — Pi), with corrections of order C(ATA^J^) or of higher order in A^~^ 

The main technical assumptions we need to prove these claims are smallness of A and 
AT. For convenience we take also the coupling 7 to the reservoirs small, as 7 = A^~^+" for 
some small a > 0. The coupling 7 is important to fix the boundary conditions for the elliptic 
equations determining T and A, but, for this, one only needs it to be bigger than A^~^. We also 
need for our analysis that uj"^ is sufficiently pinning; in fact, we will choose uj"^ = (—A + m^)^, 
with m large enough, i.e. in momentum space iu{k) = 2X]iLi(l — cosfcj) + m?. This choice 
simplifies several estimates. 

The most crucial assumption is that the space dimension has to be at least three. This is 
because of the low regularity of the collision kernels in eq. ()20|) . To prove the spectral properties 
of the linear operators we need compactness of the integral operator S in a nice enough 
space. In three dimensions this holds in a space of Holder continuous functions whereas in 
two dimensions this is not the case. Even in three dimensions, the resulting solutions have low 
regularity in momentum space which translates into long range correlations in physical space. 

We finish by comparing our equations to the standard kinetic theory (see ^2] for a discussion 
of the kinetic theory of phonon systems). In our setup, the kinetic limit is a scaling limit. 
One rescales xi to the unit interval and takes N ^ 00, X 0, while keeping R = NX'^ 
constant. Our equations then formally become, for 7 = 0, the following equation for the 
function V{x, k) = ljQ{x, k) + iJ{x, k), where x G [0, 1]: 



J"{p) = -I / dA:e-*P/2cj(p/2,A;)°sinfciJ(p/2,A;) 



(23) 



Vcj(fc) 



v^y{x,k) = RN{y) 



(24) 



with 




b{uj{k) + ijj{ki) - ijj{k'2) - uj{k2,))6{k + fci - ^2 - /C3). 
[V{ki)V{k2)V{k,) - V{k){V{ki)V{k2) + V{ki)V{k;) - V{k,)V{k,))] 



(25) 
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where the integration is over ki G [0, 2ttY. 

Equation A^(V^) = has a two parameter family of solutions corresponding to (PT|) : Vt^a{x, k) = 
uj{k)-A • Then, imposing, a posteriori, boundary conditions on equation (j211), of the form 
T(0) = Ti, T(l) = T2, A{0) = Ai, A{1) = A2, one expects, for R large, the solution of 
(|24j) to be approximately of the form Vt(x),a{x)j with T{x), A{x) solved from a scaling limit of 
our equation for them that becomes 

dx{D{T{x),A{x)){dxT{x),d^A{x)f) = 0, (26) 

where D an explicit 2x2 matrix obtained by taking the scalar product of L^^ with suitable 
vectors related to Vt(x),a{x) and orthogonal to the zero modes (see [2]). Solving with the 
prescribed boundary conditions then gives the approximate solution Vt{x),a{x) of (j^ . 

From perturbation theory one expects the full Hopf equations to reduce in the kinetic limit 
to the equation p5| . Thus our closure equations should have the same kinetic scaling limit as 
the full theory. It should be emphasized that they are more general than the kinetic equations. 
Indeed, we do not take any limit: A, r and are fixed and we prove that our solution has 
the expected properties, with precise bounds on the remainder, for A and r small and N large. 
From a mathematical point of view, the Boltzman equation has its own problems due to 
the continuum nature of the variable x. Since we work with finite, some of these unphysical 
problems are absent in our analysis. 

In conclusion, the closure equations provide an approximation to the full Hopf equations of 
the nonequilibrium state that allow us to rigorously show how the Fourier law emerges as the 
system size gets large. Moreover, this approximation should become exact in the kinetic scaling 
limit. It should provide a starting point for an eventual first principles proof of the Fourier law. 
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